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EXACT  OPERATING  CHARACTERISTICS  FOR  LINEAR  SUM  OF 
ENVELOPES  OF  NARROWBAND  GAUSSIAN  PROCESS  AND  SINEWAVE 


INTRODUCTION 

The  operating  characteristics  for  a  linear  envelope-detector  of  a 
sinewave  in  narrowband  Gaussian  noise,  followed  by  summation  of  M  independent 
envelope  samples,  were  presented  in  [1]  and  [2,  sect.  8.3].  That  approach  was 
based  upon  evaluation  of  the  first  31  moments  of  the  envelope  variate  and 
their  use  in  a  type  A  Gram-Charl ier  series  approximation,  or  in  modified 
approximations  involving  averages  over  different  numbers  of  terms  in  the 
series  [1,  pp.  758-9].  However,  there  are  possible  pitfalls  to  the  above 
approach.  First,  evaluation  of  very  low  exceedance  probabilities,  like 
10-1°,  may  be  inaccurate;  see  [1,  Fig.  1].  Second,  the  effect  of  a 
systematic  error  would  be  hard  to  detect,  if  present,  since  the  method  yields 
only  an  approximation  to  the  exceedance  distribution  function,  and  not  its 
exact  value. 

We  will  use  an  exact  approach  here,  based  upon  evaluation  of  the 
characteristic  function  of  the  envelope  detector  output,  from  which  the 
exceedance  distribution  function  can  be  precisely  evaluated  numerically 
[3,4].  In  this  fashion,  we  avoid  moment  evaluations  altogether;  we  can 
evaluate  false  alarm  probabilities  in  the  10"^  range  easily  (with  double 
precision  computer  arithmetic);  and  we  can  control  truncation  and  aliasing 
errors  to  any  desired  degree;  see  [3]  for  details.  The  results  of  [4]  can  not 
be  applied  here  because  each  independent  envelope  sample  is  the  result  of  a 
nonlinear  operation,  namely  a  sguare  root,  applied  to  a  sum  of  two  sguares  of 
Gaussian  random  variables  with  non-zero  means. 

In  the  plots  of  detection  probability  vs.  false  alarm  probability  to  be 
presented  herein,  both  abscissa  and  ordinate  use  the  same  normal  probability 
scales,  regardless  of  the  number  of  envelope  samples  M  considered.  This  allows 
for  easier  interpol ation ,  and  is  in  distinction  to  [1],  where  a  different 
false  alarm  probability  abscissa  was  used  for  each  M  [1,  pp.  759-62].  Also, 
the  parameter  employed  here  for  indexing  the  curves  is  a,  a  voltage  signal- 
to-noise  ratio  which  is  egual  to  the  ratio  of  the  sinewave  amplitude  to  the 
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rms  noise  level,  rather  than  the  dB  parameter  employed  in  [1].  This  leads  to 
curves  that  are  more  nearly  equally  spaced,  and  therefore  to  easier  and  finer 
interpolation  capability. 

Finally,  we  present  five  figures  for  the  required  input  signal-to- 
noise  ratio  per  sample  required  to  realize  specified  false  alarm  and  detection 
probabilities,  as  a  function  of  M,  the  number  of  envelope  samples  added.  The 
five  figures  correspond  to  detection  probability  P^.5,  .9,  .95,  .99,  and 
.999  respectively,  and  each  figure  contains  false  alarm  probabilities 
PpA=lO_n  for  n=l (1)8.  This  total  of  40  curves  greatly  augments  the  2 
cases  presented  in  [1,  Fig.  16]  and  [2,  Fig.  8.18]. 

A  program  for  the  evaluation  of  the  input  signal-to-noise  ratio 
required  for  a  specified  set  of  values  of  M,  PpA,  and  PQ  is  furnished, 
along  with  an  explanation  of  its  use.  In  this  fashion,  values  of  PFA’ 
and  PQ  intermediate  to  those  considered  here  can  be  easily  investigated. 
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METHOD  OF  EVALUATION 


Characteristic  Function  Details 


In  [3,4],  a  method  of  calculating  the  cumulative  and  exceedance 
distribution  functions  directly  from  a  given  characteristic  function  was 
presented-  To  utilize  those  results  here,  we  need  the  characteristic  function 
of  summation  random  variable 


x  =  e  , 
m=l 

where  e^  is  the  envelope  of  a  narrowband  filter  output  with  a  sinewave 
signal  of  amplitude  A  and  Gaussian  noise  of  power  a2.  Through  proper 
normalization,  the  probability  density  function  of  envelope  e  takes  the 
familiar  Rice  form 


;(u)  .  u  exp^  I0(cu)  for  u  >  0  , 


where  the  single  parameter 

a=7  (3) 

is  a  voltage  measure  of  signal-to-noise  ratio  per  envelope  sample.  The  power 
measure  of  signal-to-noise  ratio  per  sample  is 

2  2 
S  -  A  /2  ° 

- - ^  -  T  •  (4) 

a 

The  quantities  in  (3)  and  (4)  will  be  referred  to  as  input  signal-to-noise 
ratios,  since  they  are  per-sample  measures,  prior  to  the  summation  in  (1) 
which  yields  the  output  or  decision  variable  x. 


The  characteristic  function  corresponding  to  random  variable  e  in  (2.)  is 
given  by  Fourier  transform 


3 
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+«b  +°6  /  2  2\ 

,(  )  =  J*  du  exp(ifu)  pe(u)  =  f du  u  expM^u  -  U  ya-  )  IQ(ou) 

-a  ^  / 


and  will  be  called  the  Rice  character i Stic  function.  A  series  expansion  for 
(5)  is  developed  in  appendix  A,  and  has  been  programmed  in  double  precision 
for  numerical  use  here.  As  a  particular  special  case,  for  a=0,  no  signal,  we 
have  the  Rayleigh  probability  density  function  and  characteristic  function: 


p^0)(u)  =  u  exp(-u^/2)  for  u  >  0 


;0)(t>  -  exp(-j2/2)  j;  *  ilfl'V] 


The  latter  follows  by  use  of  [5,  3.896  3£4]  and  via  manipulation  of  the 
hypergeometric  function  series  along  with  Kumrner's  transformation  [5,  9.212  1]. 
Formula  (6)  is  particularly  attractive  numerically,  since  the  series  expansion 
of  ^  contains  all  positive  terms  except  for  one.  It  should  be  observed 

that  the  imaginary  part  of  Rayleigh  characteristic  function  f^(?)  in  (6) 
decays  very  rapidly  with?;  this  useful  feature  will  also  be  shared  by  the 
Rice  characteristic  function,  fe(?),  and  is  due  to  the  fact  that  the  odd 
part  of  the  Rice  probability  density  function  in  (2)  is  smooth  for  all  u,  and 
is  in  fact  entire  in  u,  for  any  a.  By  contrast,  the  even  part  of  the  Rice 
probability  density  function  in  (2)  has  a  discontinuous  derivative  for  real  u, 
thereby  leading  to  slow  decay  of  the  real  part  of  fe(!p. 


The  characteristic  function  of  output  variable  x  in  (1),  for 
statistically  independent  envelope  samples  [e  V  is  given  by 


fx(J)  =  [%(*>! 


in  terms  of  the  Rice  characteristic  function  (5).  This  relation  could  be  used 
directly  to  find  the  exceedance  distribution  function  of  x  according  to 
[3,  (5)— (6)3 


+«•  +«o  _ 

r  iff  fx(f) 

Qx(u)  =  J  dt  px(t)  =  2  +  J  d?  Imjexp(-iuJ)  ■—j—  . 


I 


4 
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However,  the  slow  decay  of  ^e{fx(1pi[  prompts  us  to  use  a  modified  version 
given  in  [6,  (15)]: 


+ot> 

Qx(u)  =  “  |  cos(uf)  Imjfx(f)} 


for  u  >  0 


This  form  is  applicable  to  positive  random  variables,  of  which  x,  as  given  by 
(1)  and  (2),  is  certainly  a  member. 

To  see  why  form  (9)  is  preferred  over  (8),  we  develop  (7)  as 


fx(T)  *  *  ifiirn"  ■  2.  t")  ira  [ff(r)]m  [f.ir)]""™  .  not 


where  t'r(f)  and  fi(|’)  are  the  real  and  imaginary  parts  of  Rice 
characteri Stic  function  fg(<r).  Then 


0  tv™"1  cyr>]* 


m=l 
m  odd 

contains  f.  i$)  to  at  least  the  first  power  in  all  terms,  thereby  yielding  a 
rapid  decay  with  ? . 

Development  (11)  has  been  used  to  show  why  IoM>l  decays  rapidly 
with?.  However,  when  we  employ  (9)  in  a  numerical  evaluation,  we  simply  take 
the  imaginary  part  of  the  power  in  (7),  and  do  not  use  (11)  at  all;  (11)  is  an 
alternating  series  of  large  terms  for  large  M. 


Actual  numerical  evaluation  of  (9)  proceeds  as  follows  [3]:  for  the 
Trapezoidal  rule  with  sampling  increment  A  in  ?, 


U(u) 


00 

2  1  .Cl 

—  -stm  A  +  x  “ 
tt  2  x  ^  ■  n 


cos(unA)  Im 


- 

ff  (ns)} 


I 


where  we  used  fx  0  )~1  +  iux^  as  f-*0.  Then,  restricting  the  u  values  to 
a  particular  selection. 


r  oo 

Q  (-^— — ) 

_  2 
IT 

IuxA  +  ^n  cos(2irl:,n/'0  I^{fx(nA)J 

n=l 

N-l 

=  —  ReJ>  z  exp{-i2trmn/N)  ,  (13) 

u  n=0 

where  collapsed  sequence  is  defined  as 

aO 

2o  "  ?V  *  2  7U  Im[fx'.'N4)}  • 

j  =  l 

00 

zn  =  2  nTjN  lmffx^  ^ n+J*N) A )J  for  1  <  n  <  N-l  .  (14) 

j=0 

Form  (13)  is  particularly  attractive  since  it  can  be  accomplished  via  an 
N-point  FFT.  It  can  be  shown  that  only  the  values  for  0  £  m  N/2  are  useful 
in  (13);  the  remainder  are  heavily  aliased  and  must  be  discarded.  Thus  there 
is  a  trade-off:  use  of  only  the  imaginary  part  of  f'x(J)  results  in  aliasinq 
twice  as  coarse.  However,  the  rapid  decay  of  the  imaginary  part  far  outweiqhs 
the  aliasing. 

The  summations  in  (12)  and  (14)  cannot  be  conducted  to  infinity.  Rather 
the  integral  on  f  in  (9)  is  terminated  at  limit  L,  where  the  truncation  error 
is  guaranteed  to  be  sufficiently  small.  A  trial  and  error  procedure  [3] 
yielded  the  following  rules  which  control  the  truncation  and  aliasing  errors: 

L  =  min  (9,  1 7 /VT?) , 

a  =  .12)Vfl\ 

b  =  min  (0,  -MVT727  +  V^6).  (15) 


The  inverse  VF^dependence  of  L  and  a  for  large  M  can  be  anticipated  bv 

observing  that  the  characteristic  function  ot  random  variable  x  in  (1)  then 
6 


approaches  a  Gaussian  function  with  argument  proportional  to  The  bias 

(or  shift)  b  is  added  to  random  variable  x  in  order  to  yield  a  new  random 

variable  that  remains  just  positive,  even  for  large  M;  this  allows  us  to  take 

maximum  advantage  of  the  fundamental  aliasing  interval  (0,  w/a)  in  u  in  (12) 
and  (13).  The  linear  term  (in  M)  of  b  in  (15)  is  due  to  the  mean  of  the 

Rayleigh  variate  (for  a=0)  which  is  Vtt /2v;  the  algebraic  term  inV^is  due  to 

the  fact  that  the  standard  deviation  of  random  variable  x  in  (1)  increases 
according  to  ifR*. 

In  order  to  use  this  characteristic  function  approach,  we  also  need  the 
mean  of  random  variable  x  in  (1).  Using  (2),  this  is  given  by  [5,  6.631  1] 


This  non-alternating  series  yields  accurate  values  for  the  mean. 


Special  Cases 


For  general  M,  the  characteristic  function  approach  described  above  must 
be  used.  However,  for  M  =  1  and  2,  closed  form  expressions  for  the  false 
alarm  and  detection  probabilities  are  possible.  Specifically,  from  (1)  and 
(2) ,  for  u  >  0, 


And  for  M  =  2,  the  false  alarm  probability  can  be  determined  by  convolving  two 
Rayleigh  probability  density  functions  of  the  form  of  (6),  to  give,  for  u  _>  0, 
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Pc  a  =  exp(-u  y2)  +y?u  exp(-u^/4) 


[%) '  * 


for  M  =2.  (18 


Here,  $  is  the  cumulative  distribution  function  of  a  normalized  Gaussian 


random  variable: 


f(u)  =  J  dt  (2ir)-1^2  exp(— 12/2 ) 


The  detection  probability  of  random  variable  x  in  (1)  is  not  available  in 
closed  form  for  M  >  1. 

Asymptotic  Performance  for  Large  M 

For  large  M,  decision  variable  x  in  (1)  is  approximately  Gaussian.  The 
mean  of  x  was  given  in  (16);  a  similar  approach  for  the  mean  square  of  x 
yields  the  variance  as 


a2  =  M a2  =  M(2+a2-„2)  . 


The  probability  density  function  of  x  is  then  approximately 


with  exceedance  distribution  function 


qx(u)  s  ) = if— =— 

\°x/ 


For  input  signal-to-noise  ratio  S/N=*0,  we  have  a=0  from  (4),  and  (22), 
(16),  and  (20)  specialize  to 
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On  the  other  hand,  for  S/N>0,  (22)  yields  the  detection  probability  P^.  We 
now  use  the  inverse  function  |>  to  definition  (19)  and  solve  (23)  and  (22) 
according  to 


MAS72-U  f,  ,  M“e'u  , 
— 7— —  -  SO’fa'  •  7£T~  -$(PD> 


Eliminating  threshold  u  in  (24),  we  have 


»e*p0>  ■ 


But  also,  for  large  M,  the  required  per-sample  input  signal-to-noise  ratio  « 
will  be  small,  giving 


~]/T ifi(  2* i;  ■  z)  ay?(i+4 J  * 


2  ~  .  2  2  0  n  .  , 

ae=2a_ue2*2  '  (26) 

Substituting  these  results  in  (25)  and  solving  for  a,  we  have  the  required 
per-sample  input  signal-to-noise  ratio  measures  for  large  M  in  the  alternative 
forms 


/,  \l>*  a1'2  _  1/2 

77T'  1-446  ^p-  • 


• 


>,  1.045  -jr. 


dB  =  10  log  |  a  10  log(2p^J+10  logU)-5  log(M)  =  .193+10  log(s)-5  loq(M),  (27) 


where  the  single  parameter 


8  =  $(Pn)  -  $(PF  J 


9 


RESULTS 


For  a  given  value  of  M,  the  output  variable  in  (1), 

*  -  •  (29> 
m=l 

will  exceed  threshold  u  with  false  alarm  probability  PpA  when  signal-to-noise 
ratio  «  is  zero.  That  is 

PpA  =  Prob(x>u  l  c=0;  M).  (30) 

For  specified  values  of  M  and  PpA,  this  relation  can  be  solved  numerically 
for  u;  the  values  of  normalized  threshold  u/M  are  listed  in  table  1  for 
M=2n,  n=0(l)13  and  for  PFA=10'n,  n=l(l)8. 

The  detection  probability  depends  on  threshold  u,  M,  and  signal-to-noise 
ratio  a(>0): 

Pq  =  Prob(x>u  |  a;  M).  (31) 

For  specified  values  of  M,  P^,  and  u,  this  relation  can  be  solved 
numerically  for  the  required  input  signal-to-noise  ratio  a.  When  the 
threshold  results  in  Table  1  are  employed,  the  results  yield  the  required 
input  signal-to-noise  ratio  for  specified  false  alarm  probability  and  detection 
probability  at  a  particular  M.  These  are  plotted  in  figures  1-5  for 

PD  =  .5,  .9,  .95,  .99,  .999,  (32) 

respectively.  The  abscissa  is  loq^M,  and  the  ordinate  is  in  decibels,  as 
defined  in  (27).  The  fit  of  (27)  is  very  qood  for  large  M,  especially  for 
the  larger  PpA  values.  These  results  in  figures  1-5  greatly  extend  the  one 
in  [1>  Fig.  16]  and  [2,  Fig.  8.18]. 
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Table  1.  Normalized  Thresholds  Required  for  Specified  M  and  PpA 


\p 

M\FA 

IE-1 

IE-2 

IE-3 

1 E  -  4 

1 

2.14596603 

3.03485426 

3,71692219 

4,29193205 

2 

1.87154046 

2.46578168 

2,92459903 

3 . 3  1  3  7  2 7  9 

4 

1.68491649 

2 . 08494224 

2 , 39281962 

2.65432267 

8 

1 ,55592564 

1 .82779134 

2.03544098 

2,21134522 

16 

1 ,46605729 

1,65246898 

1 . 79362769 

1 ,91266565 

32 

1,40314416 

1  .53192213 

1 ,62866385 

1 .70984877 

64 

1.35896377 

1 ,44846093 

1 .51524477 

1,57104117 

128 

1.32787317 

1 ,39035933 

1 ,43673968 

1 ,47534630 

256 

1,30596258 

1,34974198 

1.38210498 

1 .40896493 

512 

1.29050601 

1 .32125803 

1 ,34392160 

1 ,36268942 

1024 

1.27959472 

1.30123656 

1 .31715039 

1 .33030674 

2048 

1.27188832 

1 ,28713956 

1.29833595 

1 ,30758095 

4096 

1 .26644357 

1 .27720181 

1.28509047 

1 ,29159844 

8192 

1.26259580 

1.27018998 

1 .27575385 

1 .28034098 

FA 


IE-5 


IE-6 


IE-7 


IE-8 


1 

4.79852591 

5.25652177 

2 

3.65817649 

3.97074674 

4 

2.88639585 

3,09755766 

8 

2,36734857 

2.50933650 

16 

2,01795589 

2.11363367 

32 

1,78141625 

1,84629005 

64 

1 ,62006566 

1 ,66438962 

128 

1,50917003 

1.53967893 

256 

1 .43244246 

1  .  45357776 

512 

■  1.37906412 

1,39378248 

1024 

1 .34176981 

1 ,35206123 

2048 

1.31562790 

1.32284604 

4096 

1.29725887 

1,30233301 

8192 

1.28432858 

1,28790149 

5  *  67769243 
4  *  25904998 
3.29282208 
2 . 64073862 
2.20209577 
1  .90615996 
1 .70520835 
1.56771937 
1 .47297026 
1.40726893 
1 .36148146 
1 . 32944798 
1 .30697131 
1.29116614 


6.06970852 
4.52806135 
3.47544423 
2.76376208 
2.28487698 
1  .96210527 
1 .74328423 
1 .59383081 
1 .49100181 
.1 .41979378 
1 ,37022180 
1 ,33556910 
1 ,31126956 
1 ,29419029 


Figure  1.  Required  Input  S/N  for 


Figure  5.  Required  Input  S/N  for 
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The  results  in  figures  1  through  5  only  cover  a  selected  set  of  detection 
and  false  alarm  probability  values.  A  more  complete  description  is  afforded 
by  the  receiver  operating  characteristics,  namely  detection  probability  vs. 
false  alarm  probability,  with  signal-to-noise  ratio  as  a  parameter.  In 
figures  6  through  19  are  given  these  operating  characteristics  for 

M  =  1,  2,  4,  8,  16,  32,  64,  128,  256,  512,  1024,  2048,  4096,  8192,  (33) 

respectively.  The  false  alarm  probability  covers  the  range  10“^  to  .5, 
while  the  detection  probability  covers  10“^  to  .999.  Both  abscissa  and 
ordinate  in  these  figures  employ  the  inverse  function  to  the  Gaussian 
cumulative  distribution  function  5  defined  in  (19);  thus,  a  truly  Gaussian 
random  variable  would  plot  as  a  series  of  equally  spaced  parallel  straight 
lines  (with  parameter  a).  Observe  that  the  curves  are  nearly  equally  spaced 
with  parameter  a,  except  for  very  small  a,  where  the  nonlinear  envelone 
operation  causes  small  signal  suppression  and  a  crowding  together  of  the 
curves. 

If  the  decision  variable  x  is  presumed  Gaussian,  and  the  ooerating 
characteristics  overlayed  on  the  exact  results  in  figures  6-19,  it  is  found 
that  the  two  sets  of  curves  for  M=8192  are  virtually  identical  in  the  range  of 
ppA  and  Pq  plotted.  However,  for  M=16,  the  Gaussian  approximation  is 
somewhat  optimistic;  for  example,  the  exact  curve  for  «=2.75  is  well- 
approximated  by  the  Gaussian  approach  for  a=2.62.  For  small  M,  the  Gaussian 
approximation  is  overly  optimistic  for  small  PpA;  however,  the  two  sets 
cross  near  Pp^=.5,  which  is  not  a  practical  range  of  interest  anyway. 


IS 


PROBABILITY  OF  DETECTION  PD 


PROBABILITY  OF  DETECTION  PD 
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Figure  13.  Receiver  Operating  Characteristics  for  M=128 
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Figure  15.  Receiver  Operating  Characteristics  for  M=512 


IE-10  1 E - 7  IE-4  1 E -  3  IE-2  IE-1 
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Figure  19.  Receiver  Operating  Characteristics  for  M=8192 
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SUMMARY 

A  method  for  exact  evaluation  of  the  exceedance  distribution  function,  of 
a  linear  sum  of  M  envelopes  of  a  narrowband  Gaussian  process  and  sinewave,  has 
been  utilized  to  determine  the  receiver  operating  characteristics  for  a  wide 
range  of  values  of  M  and  signal-to-noise  ratio.  Also,  the  required  input 
signal-to-noi se  ratio  vs.  M  has  been  determined  for  a  selected  set  of  false 
alarm  and  detection  probabilities.  Programs  are  also  supplied  by  which  other 
values  of  the  various  parameters  can  be  investigated  by  the  user. 

Agreement  between  the  current  results  and  those  in  [1,2]  is  very  good 
over  the  range  of  common  values  plotted.  For  M  larger  than  8192,  the 
approximation  given  in  (27)  and  (28)  is  recommended,  since  the  summation 
variable  is  then  well  represented  by  a  Gaussian  random  variable. 
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APPENDIX  A.  DERIVATION  OF  RICE  CHARACTERISTIC  FUNCTION 

The  normalized  probability  density  function  of  a  Rice  random  variable  was 
given  in  (2)  as 


Pe(u)  =  u  exp 


f  ^ .. 


(au)  for  u  >  0 


( A— 1 ) 


The  corresponding  characteristic  function  is 


+  7(2  2\ 

fe(?)  =  J*du  exp(ifu)  pe(u)  =  du  u  exp  ^iju  -  -u  IQ(au)  = 

=  exp(-r)  5  I—  [* du  u2n+1  exp(  ifu-u2/2)  , 

„  (n!  j  J 


(A-2) 


where  we  have  expanded  Iq  in  a  power  series  [5,  8.447  1]  and  defined  power 
signa 1-to-noise  ratio 


r  =  ac/2  . 


( A— 3 ) 


(If  desired,  a  power  series  in  f  could  be  developed  by  expanding  exp(ifu)  in  a 
power  series  instead  of  IQ.) 

We  define 


C  (?)  =  — 2  fdu  u2n+1  exP^Tu  - 

0  ^(n!)2  \ 

and  get  the  characteristic  function  series 

fe(f )  =  exp(-r)  rn  C^f) 
n=0 


In  order  to  get  a  recurrence  on  Cn(f),  we  also  define 


12)  for  n  >  0  ,  (A-4) 


(A— 5 ) 


•O 

.(f)  =  j"dw  wk  expf.ijw  -  w2/2)  for  k  >  0  , 


(A— 6 ) 
A-l 


‘4 

n 


•* .  •*. 

.  •  m.  *  _  m  m 
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for  then 


c„<r>  - 


2n{n :)2 


By  integrating  by  parts  on  (A-6),  there  follows 


8k  =  ^Sk-1  +  ( k-1 ) 8^-2  for  k  >  1 


This  recurrence  can  be  started  with  [5,  3.896  3^4] 


i0  =e,p(-f2/2)[|f.  ,? 


By  looking  at  three  adjacent  terms  of  recurrence  (A-8),  we  can  generate 
the  alternative  recurrence 


B.  =  (2k-3-/)B,2  -  ( k— 2 )  (k-3)B._. 


By  means  of  (A-7),  this  translates  into 


for  n  >  2  .  (A- 


Starting  values  are  (via  manipulation  of  hypergeometric  series  and  Kummer's 
transformation)  expressable  as 

C0  =  exP(-^/2)  |^1F  1  ^  2;  2’  f“)  +  » 

Cj  =  exp(-?2/2)  LFj(-  f;  pj  *  i|f  (3-{2)-f].  (A- 

Each  of  the  series  for  consists  of  terms  of  the  same  polarity,  except 
for  one  term,  and  are  therefore  useful  for  obtaining  very  accurate  initial 
values.  Cq  is  the  characteristic  function  of  the  Rayleigh  probability 
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density  function.  Relations  ( A— 1 1 ) — ( A— 1 2 )  constitute  recurrences  on  both  the 
real  and  imaginary  parts  of  Cn. 

It  was  found  that  the  terms  exp(-r)  r11  in  (A-5)  became  very  large  for 
large  n,  while  the  Cn(^)  terms  became  very  small.  In  order  to  avoid 
overflow  and  underflow,  we  defined  the  total  term 

An  =  exp(-r)  rn  Cn(f)  .  (A-13) 

Reference  to  ( A— 1 1 )  readily  yields  the  recurrence  on  A  ,  and  (A-12) 
furnishes  corresponding  obvious  starting  values  for  Aq  and  A^. 


A-3/A-4 
Reverse  Blank 


TR  7117 


APPENDIX  B.  DESCRIPTION  OF  PROGRAMS  AND  LISTINGS 


Overview 


Information  obtained  via  evaluation  of  the  Rice  characteristic  function 
may  be  displayed  in  three  formats. 

FORMAT  1:  Display  PD  vs.  PFA 

The  user  defines  the  number  of  samples  M  and  the  range  of  values  for 
alpha,  a  voltage  signal-to-noise  ratio  measure.  An  algorithm  then  utilizes 
the  Rice  characteristic  function  for  alpiia=0  and  for  the  alphas  specified  by 
the  user.  This  results  in  the  production  of  a  threshold  vs.  PFA  and  M 
(alpha=0)  and  threshold  vs.  PD  and  M  (alpha>0)  tables.  These  two  tables  are 
stored  on  an  output  file.  For  each  user-defined  M,  a  plot  routine  displays  PD 
vs.  PFA  for  the  set  of  user-defined  alphas. 

FORMAT  2:  Display  SNR  vs  M 

The  user  supplies  the  input  which  specifies  a  PD.  The  algorithm  then 
solves  for  the  threshold  values  corresponding  to  PFA=10**(-IPFA) , 

(IPFA=1, . . ,8)  and  M=2**IM,  ( IM=0, . . . ,13)  and  alpha=0.  A  root  finding 
technique  is  then  employed  to  solve  for  the  SNR  defined  by  a  threshold  value 
and  user-defined  PD.  An  SNR  is  found  for  each  threshold  value.  The  results 
are  stored  in  an  output  file.  A  plot  routine  displays  the  required  SNR  vs.  M 
for  PFA=10**(-IPFA) ,  (IPFA=1,2, . . . ,8) . 

FORMAT  3:  Print  SNR 

The  user  specifies  a  value  for  PD,  PFA,  M.  The  program  solves  for  the 
threshold  corresponding  to  PFA  and  M.  A  root  finding  technique  is  then 
employed  to  determine  the  SNR  corresponding  to  this  threshold  and  user-defined 
PD  and  M.  The  results  are  printed. 
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Description  of  Input 


Inputs  to  the  program  consist  of  cards  which  either  specify  values 
(PARAMETER  CARDS),  activate  the  reading  of  tabularized  values  (TABLES),  assiqn 
files  (FILE  NAME  CARDS),  process  data  (COMMAND  CARDS),  or  specify  a  plot 
device  (PLOT  DEVICE  CARDS).  The  basic  format  of  a  card  is 

CARD  NAME  =  value  units 

where  CARD  NAME  is  an  alphanumeric  expression  from  Tables  2-6.  The 
alphanumeric  must  begin  in  column  1,  value  is  a  floating  point  or  integer 
number,  and  units  is  an  alphanumeric. 

Parameter  cards,  file  names,  and  tables  constitute  the  data  upon  which 
commands  operate.  If  two  cards  with  the  same  name  specify  different  data, 
then  the  last  entry  overrides  the  other. 

For  the  programmers  convenience,  FORTRAN  variable  names  associated  with 
file  names  or  parameters  may  be  located  in  the  Tables  2  through  6.  Since 
input  and  values  stored  represent  the  same  physical  quantity,  it  is  convenient 
to  refer  to  both  in  this  paper  by  the  same  variable  name.  The  convention 
adopted  is  to  express  the  variable  by  the  lower  case  letters  and  reserve  upper 
case  letters  for  constants. 

Parameter  Cards 


Parameter  cards  are  used  to  specify  an  axis  length  or  assiqn  a  range  of 
values  to  a  parameter.  These  cards  are  shown  in  Table  2.  For  example, 

NUMBER  OF  SAMPLES  MINIMUM  =  1. 

NUMBER  OF  SAMPLES  MAXIMUM  =  8192. 

NUMBER  OF  SAMPLES  FACTOR  =  2. 


implies  that  the  program  will  process  data  for  M=l,2,4,8,16, . . . . ,4096,8192. 
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Table  Cards 

A  table  card  contains  the  values  that  are  to  be  assigned  to  a  variable. 
The  last  card  that  must  appear  in  a  table  is  an  EOF  card.  This  card 
terminates  the  reading  of  the  table.  Table  cards  exist  for  PD  and  PFA  only. 

A  list  of  the  table  cards  appears  in  Table  4.  For  example, 

PROBABILITY  OF  DETECTION  TABLE 

.5 

.7 

.99 

EOF 

This  table  assigns  values  of  .5,  .7,  .99  to  PD. 

Files  Cards 

A  file  card  allows  for  dynamic  assignment  of  all  mass  storage  files. 

This  is  accomplished  by  linking  internal  FORTRAN  unit  numbers  to  files  during 
execution.  The  file  card  is  shown  in  Table  4.  Two  of  the  three  algorithms 
use  files.  They  are 

Display  PD  vs  PFA  :  A  file  is  used  to  store  output. 

Display  SNR  vs  M  :  A  file  is  used  to  store  output. 

For  example, 

OUTPUT  FILE  =  PDFILE 

directs  the  output  of  a  program  to  a  file  called  PDFILE. 

Command  Cards 

Command  cards  are  used  to  compute,  plot,  or  terminate  a  run  stream. 

Command  cards  are  given  in  Table  5. 
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Plot  Device  Cards 

Plot  device  cards  direct  the  plot  output  to  either  a  TEKTRONIX,  FR80,  or 
a  CALCOMP  plotter.  The  cards  necessary  for  that  operation  are  shown  in 
Table  6. 

Examples  of  Output 

Example  1:  Display  PD  vs  PFA 

The  input  deck  for  the  first  example  appears  in  Table  7.  This  deck 
designates  that  PD  vs.  PFA  data  will  be  computed  for  M=1  and 
alphas. 5, 1.0, 1.5, . . . ,9.5.  The  output  is  stored  on  a  file  called  FILE1.  The 
plot  corresponding  to  the  data  is  shown  in  figure  6.  The  second  half  of  the 
run  stream  computes  PD  vs.  PFA  data  for  M=2  and  alpha=0. , .4, .8, . . . ,7.2.  The 
output  is  stored  in  FILE2.  The  plot  of  the  data  appears  in  figure  7. 

Example  2:  Display  SNR  vs.  M 

The  input  deck  for  the  second  example  appears  in  Table  8.  The  first  half 
of  the  input  deck  designates  that  the  SNR  vs.  M  plots  will  be  computed  for  a 
value  PD=.5.  The  output  is  displayed  in  figure  1.  The  parameter  cards 
specify  that  the  axis  will  be  scaled  as  follows:  -19  DB  (minimum),  13  DB 
(maximum),  2  DB  (increment),  and  5  inches  long  for  the  SNR  axis  and  6.86 
inches  long  for  the  number  of  samples  axis.  It  should  be  noted  that  the 
limits  for  the  number  of  samples  axis  are  predefined  by  the  program  to  be  1 
(minimum),  8192  (maximum),  2  (factor).  The  output  is  stored  in  a  file  called 
PDFIll.  The  second  half  of  the  run  stream  computes  SNR  vs.  M  for  a  value 
PD=.9.  The  axis  limits  for  SNR  were  changed  to  -17  DB  (minimum),  15  DB 
(maximum),  2  DB  (increment).  Alpha  curves  were  computed  for 
alpha=0. , .4, .8, . . . ,7.2.  This  output  is  stored  in  file  PDFIL2.  A  plot  of  this 
data  appears  in  figure  2. 

Example  3:  Print  SNR 

The  input  deck  for  the  third  example  appears  in  Table  9.  The  output 
appears  in  Table  10. 
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TABL  E  2  ►  PARAMETER  CARDS 


INPUT  CARDS  UNITS 


SNR  AXIS  LENGTH  *  snraxs  IN 

SAMPLE  AXIS  LENGTH  »  sinpaxs  IN 

PD  AXIS  LENGTH  =  pdaxs  IN 

PFA  AXIS  LENGTH  *  pfssxs  IN 

SNR  MINIMUM  =  snritin  DB 

SNR  MAXIMUM  =  snrmax  DB 

SNR  INCREMENT  -  snrinc  DB 

ALPHA  MINIMUM  «  alpmin 
ALPHA  MAXIMUM  =  alpmax 
ALPHA  INCREMENT  *  alpine 


NUMBER  OF  SAMPLES  MINIMUM  smpmin 
NUMBER  OF  SAMPLES  MAXIMUM  -  smpmax 
NUMBER  OF  SAMPLES  FACTOR  ~  smpf ct 


TABLE  3 »  TABLE  CARDS 


INPUT  CARDS  VARIABLE 


PROBABILITY  OF  DETECTION  TABLE  PD 

PROBABILITY  OF  FALSE  ALARM  TABLE  PFA 


TABLE  A,  FILE  CARDS 


INPUT  CARDS 


OUTPUT  FILE  ~  name 


TABLE  5,  COMMAND  CARDS 


INPUT  CARDS 
RUN  MAIN 

COMPUTE  PD  VS  PFA 
COMPUTE  SNR  VS  M 
PLOT  PD  VS  PFA 
PLOT  SNR  VS  M 
END 
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TABLE  6.  PLOT  DEVICE  CARDS 


INPUT  CARDS  OPTIONS 


BAUD  RATE  =  960. 

PLOT  DEVICE  =  device  FR80 » TEKTRO > CALCOMP 

RESET  PLOT  DEVICE 


TABLE  7,  SAMPLE  INPUT  DECK  FOR  PD  VS  PFA 

RUN  MAIN 

BAUD  RATE  =  960. 

PLOT  DEVICE  =  TEKTRO 

RESET  PLOT  DEVICE 

PD  AXIS  LENGTH  »  6.86  IN 

PFA  AXIS  LENGTH  =  5.  IN 

OUTPUT  FILE  -  FILE1 

NUMBER  OF  SAMPLES  MINIMUM  1 

ALPHA  MINIMUM  =  .5 

ALPHA  MAXIMUM  »  9.5 

ALPHA  INCREMENT  =  .5 

COMPUTE  PD  VS  PFA 

PLOT  PD  VS  PFA 

OUTPUT  FILE  =  FTLE2 

NUMBER  OF  SAMPLES  MINIMUM  =  2 

ALPHA  MINIMUM  -  0. 

ALPHA  MAXIMUM  *  7.2 
ALPHA  INCREMENT  ■  ,4 
COMPUTE  PD  VS  PFA 
PLOT  PD  VS  PFA 
END 
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TABLE  8 .  SAMPLE  INPUT  DECK  FOR  SNR  vs  M 

RUN  MAIN 

BAUD  RATE  =  960. 

TEMPORARY  FILE  =  FALSE 

PLOT  DEVICE  »  TEKTRO 

RESET  PLOT  DEVICE 

OUTPUT  FILE  =  PDFIL1 

SNR  MINIMUM  =  -19.  DB 

SNR  MAXIMUM  =  13.  DB 

SNR  INCREMENT  =  2.  DB 

SNR  AXIS  LENGTH  =  5.  IN 

SAMPLE  AXIS  LENGTH  -  6.86  IN 

PROBABILITY  OF  DETECTION  TABLE 

.5 

EOF 

COMPUTE  SNR  VS  M 

PLOT  SNR  VS  M 

OUTPUT  FILE  *  PDFIL2 

SNR  MINIMUM  =  -17.  DB 

SNR  MAXIMUM  «  15.  DB 

SNR  INCREMENT  *  2.  DB 

PROBABILITY  OF  DETECTION  TABLE 

.9 

EOF 

COMPUTE  SNR  VS  M 
PLOT  SNR  VS  M 
END 


TABLE  9.  SAMPLE  INPUT  DECK  FOR  PRINTING  SNR 


RUN  MAIN 

PROBABILITY  OF  DETECTION  TABLE 

.5 

.9 

EOF 

PROBABILITY  OF  FALSE  ALARM  TABLE 
.  1 

.001 

EOF 

NUMBER  OF  SAMPLES  MINIMUM  =  1. 
NUMBER  OF  SAMPLES  MAXIMUM  =  2048. 
NUMBER  OF  SAMPLES  FACTOR  =•  2. 
PRINT  SNR 
END 
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This  section  contains  a  listing  of  three  master  programs  and  associated 
subroutines.  Subroutines  which  read  input  and  plot  the  output  have  been 
omitted.  Table  11  contains  a  list  of  the  subroutine  names  and  a  brief 
description  of  the  pertinent  subroutines. 
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TABLE  11.  DESCRIPTION  OF  SUBROUTINES 


NAME 


CMPDVA 

CMPSVS 

PRTSNR 

FFT 

RDC 

FNPD 

FNPF 

FNF11 

RICE 

FNIPHI 

DIST 


DESCRIPTION 


MASTER  PROGRAM  FOR  COMPUTING  PD  US  PFA 
MASTER  PROGRAM  FOR  COMPUTING  SNR  VS  M 
MASTER  PROGRAM  FOR  COMPUTING  AND  PRINTING  SNR 
COMPUTES  THE  FAST  FOURIER  TRANSFORM  OF  A  FUNCTION 
COMPUTES  AN  APPROXIMATE  S/N  FOR  A  GIVEN  PD>  PFAs-  M 
< SEE  REF  7) 

COMPUTES  THE  PROBABILITY  OF  DETECTION  FOR  A  GIVEN 
M,  S/N ?  AND  THRESHOLD 

COMPUTES  THE  PROBABILITY  OF  FALSE  ALARM  FOR  A  GIVEN 
M  AND  THRESHOLD 

COMPUTES  THE  CONFLUENT  HYPERGEOMETRIC  FUNCTION 
COMPUTES  THE  CHARACTERISTIC  FUNCTION  OF  A  RICE 
VARIATE 

COMPUTES  THE  INVERSE  OF  THE  CUMULATIVE  GAUSSIAN 
DISTRIBUTION 

COMPUTES  THE  EXCEEDANCE  DISTRIBUTION  FUNCTION  FOR 
A  GIVEN  M  AND  S/N 


1 
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SUBROUTINE  FFT (  N  *  X  >  Y  ) 

IMPLICIT  DOUBLE  PRECISION  <A-H»0~Z> 

DIMENSION  C(0J256)  *X(0*1023)  .<Y(051023)  >  L  <  0  *  9 ) 
DATA  PI/3* 1415926535897932400/ 

T=2.D0*PI/N 

Jl=N/4 

DO  100  J  =  0,J.t 
C< J)=DCOS<T*DFLOTJ< J) > 

100  CONTINUE 

Nl=N/4 
N2=N1 +1 
N3=N2+ 1 
N  4  =  N  3  +  N 1 

L2=JIDINT(1 , 4427D0*DL0G<DFLDTJ<N>)+.5D0> 

DO  600  1 1  =  1  »  L2 
I2=2*#<L2-I1 > 

I3=2D0*I2 
14  =N/I3 

DO  500  15=1 » 12 

16=14# ( 15-1 ) +1 

IF (  1 6  *  LE  *  N2  )  60  TO  350 

V6  =  -C ( N4-1 6-1 ) 

V7=-C<I6-N1-1 > 

GO  TO  375 
350  V  6  =  C  ( 16-1 ) 

U7  =  -C  ( N3-I6-.1 ) 

375  L.3  =  N- 1 3 

DO  400  I  7  =  0  »  L3  > 13 

18=17+15 

19=18+12 

V8=X( 18-1 )-X< 19-1 ) 

V9  =  Y< 18-1 >-Y< 19-1 ) 

X< 18-1 >=X< 18-1 )+X< 19-1  ) 

Y( 18-1 )=Y( 18-1 >+Y( 19-1 > 

X<  19-1  )=V6*V  8~U7*V9 
Y< 19-1 ) - V  6  #  V  9  +  V  7  # V  B 
400  CONTINUE 
500  CONTINUE 
600  CONTINUE 


1 1  =  L  2  + 1 

DO  700  12  =  1 r  1  0 
L ( 12-1 )=1 ►DO 

I F (  I2.GT.L2  )  GO  TO  700 
L< 12-1 >=2**< 11-12) 

700  CONTINUE 
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ICO  =  L ( 0 ) 
IC1=L<  1 ) 
IC2=L<2) 
IC3--L.  ( 3 ) 
I C  4  =  L  (  4  ) 
I C5  =  L ( 5  > 
IC6=L(6> 
IC7  =  L ( 7 ) 
IC8  =  L  <  8 ) 
IC9=L<9> 


K-l 

DO  1900  I 1=1 » IC9 
DO  1800  12=11 , IC8, IC9 
DO  1700  13=12, IC7, IC8 
DO  1600  1 4  =  1 3  * IC6 , I C7 
DO  1500  15=14, IC5, IC6 
DO  1400  16=15, IC4  > IC5 
DO  1300  17=16, IC3, IC4 
DO  1200  18=17 » IC2  ,  IC3 
DO  1100  1 9= 1 8 , 1 C 1 » IC2 
DO  1000  110=19, ICO, IC1 
J=I10 

IF <  K.GT.J  )  GO  TO  900 
A=X ( K- 1 ) 

X  (  K-l ) =X  <  J-l ) 

X< J-l )=A 
A=Y ( K-l > 

Y<K-1 >=Y< J-l ) 

Y ( J-l ) =A 
900  K=K+1 

1000  CONTINUE 
1100  CONTINUE 
1200  CONTINUE 
1300  CONTINUE 
1400  CONTINUE 
1500  CONTINUE 
1600  CONTINUE 
1700  CONTINUE 
1800  CONTINUE 
1900  CONTINUE 


RETURN 
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SUBROUTINE  FNPD  ( ALPHA » V  *  AM  »  AL  » AD » ABS>  PD ) 
IMPLICIT  DOUBLE  PRECISION  <A-H»0-Z> 

DATA  PI/3. 141 59265358979324D0/ 

SNR= ,5D0*ALPHA*ALPHA 

CALL  FNF11 ( 1 .  5D0  > 1 .  DO  >  SNRuFl 1 ) 

FAC=DSQRT ( . 5D0*PI ) *EXP < -SNR ) *F 1 1 
AMUY=AM*FAC+ABS 
AM2-  AM/2 .  DO 
VD  =  V*AD 

EXC=.5*AD*AMUY 
NS1~JIDINT (AL/AD) 

DO  100  NS  =  1 >NS1 
XI=AD*NS 

CALL  RICE  <  XI >SNRrFRxFI ) 

A=DATAN2(FI»FR> 

FYI^  DSIN<AM*A  +  ABS*XI)*(FR*FR  +  FI*FI)**AM2 
ADD=FYI*DCOS< VD*DFLOTJ < NS > >/DFLOTJ(NS> 
EXC=EXC+ADD 
100  CONTINUE 

PD=2.D0*EXC/PI 

RETURN 

END 


SUBROUTINE  FNPF <  V  »  AM » AL  *  AD  * ABS»  PF ) 
IMPLICIT  DOUBLE  PRECISION  <A-H»0-Z> 
DATA  PI/3. 141592A535B979324D0/ 


FAC=DSQRT ( ,5D0*PI ) 

AMUY=AM*FAC+ABS 
AM2=AM/2 . DO 
VD=V*AD 

EXC;:  .5*AD*AMUY 
NS1"JIDINT (AL/AD) 

DO  100  NS-^IjNSI 
XI ~AD*NS 
X2=.5D0*XI*XI 
E=EXP ( -X2 ) 

CALL  FNF11  < - . 5D0 »  ,5D0)»X2^F11  ) 

FR=E*F 1 1 
FI~E*FAC#XI 
A~DATAN2  <  F I >FR) 

FYI^DSIN( AM*A+ABS*XI ) * < FR*FR+FI *FI >**AM2 

ADD  =  FYI*DCOS<  VD*DFLOT.J<NS>  )  /DFLQT  J  <  NS ) 

EXC=EXC+ADD 

CONTINUE 

PF=2.D0*EXC/FI 


RETURN 

END 


SUBROUTINE  RDC  ( AM  >  PF  >  PD  >  ALPHA ) 

IMPLICIT  DOUBLE  PRECISION  <A-H>0-Z) 

A=DLOG< . 62D0/PF ) 

B=DLOG  <  PD/ ( 1 .DO-PD)  > 

FACT ”6  » 2D0  +  4 . 54D0/DS0RT < AM  + . 4 4DO ) 

SNRDB=-5.D0*DL0G10<AM>  +  DLOG 1 0  <  A  +  .  1 2D0*A*B+ 1 . 7D0#B  )  *FACT 
ALPHA=DSGRT(2.D0*10 ,DO**< ► 1 DO#SNRDB )  ) 

RETURN 

END 


SUBROUTINE  FNF 1 1 < A , B , X  * F 1 1  ) 

IMPLICIT  DOUBLE  PRECISION  (A-HrO-Z) 

F 1 1  ~  1  .DO 
T~1 » DO 

DO  100  K  ■=  1  x  300 
U=K~1 

T=T*( A+U)*X/< ( B+U ) *K ) 

FI 1 -F 1 1 +T 

IF (  DABS <T) fLE»DABS(Fll >  *1 >0-18  )  GO  TO  200 

100  CONTINUE 
PRINT  101 

101  FORMAT (2X» '300  TERMS  IN  FNF1 1 ' > 

200  CONTINUE 

RETURN 

END 


SUBROUTINE  FNIPHI <  X » PHI ) 

IMPLICIT  DOUBLE  PRECISION  (A-Hi-O-Z) 

100  Y=DMAX1 <X» 1 ♦D-12) 

Y=DMIN1<Y»1 ,D0-1 »D-12 ) 

D=X- *  5D0 

IF (  DABS ( D ) » GT .  .OIDO  )  GO  TO  250 

PHI=2 ► 506628274 63D0#D*  < 1  *  D0+D#D#1 *  0471 9755120D0 ) 

GO  TO  300 
250  PHI=Y 

IF  <  Y . GT .  . 5D0  )  PHI=»5D0-(Y-*5D0) 

PHI-DSGRT <-2*D0#DL0G<PHI )  ) 

T=1 * DO+PHI # < 1 ,43278BD0+PHI#< . 189269D0+PHI* ♦ 001 308D0 >  ) 
PHI=PHI-(2.515517D0+PHI*( ♦  802853D0+PHI* *  01 0328D0 > )/T 
IF (  Y  »LT .  • 5D0  )  PHI=-PHI 
300  RETURN 
END 
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SUBROUTINE  DIST < AH » ALPHA t HF t X » Y ) 
IMPLICIT  DOUBLE  PRECISION  <A-H»0-Z> 
DIMENSION  X<0‘*1023>,  Y<0tl0?3> 
COMMON  /PDVPF/AL>ADi-ABS 
DATA  PI/3* 141 59 265358979324 DO/ 


SNR=  ►  5D0*ALPHA#ALPHA 

CALL  FNF1 1 < 1 *5D0fl *DO»SNR»Fll ) 

AMU=DSGRT<  .  5D0*P  I )  *DEXP  ( --SNR  >  *F1 1 

AMUS=AM*AMU+ABS 

AM2=AM/2  *  DO 

DO  100  I-Or 1 023 
X  ( I )  =0  *  DO 
Y  ( I )  =0  *  DO 
100  CONTINUE 

X<0>=.5D0*AMUS#AD 
NS1=JIDINT  <  AL/AD ) 

DO  1000  NS=1 > NS1 
XI=AD*NS 

CALL  RICE(XI »SNR»U»V) 

T=DATAN2(V>U) 

FI=DSIN<AM*T+ABS*XI)*<U*U+V*V)**AM2 

MS=JMOD(NS»MF) 

X(MS)=X(MS)+FI/NS 
1000  CONTINUE 

CALL  FFT ( MF  » X  >  Y ) 

FAC- 2 . DO/PI 
KS1 =MF/2  *  DO 
DO  2000  KS~0»KS1 
T=X(KS>*FAC 
X ( KS ) =1 .DO-T 
Y<KS)=T 

2000  CONTINUE 

RETURN 

END 
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SUBROUTINE  RICE  <  X* SNR  *  FR  *  FI ) 
IMPLICIT  DOUBLE  PRECISION  <A-H»0-Z) 
DATA  PI/3. 14159265358979324D0/ 


X2=.5D0*X*X 

Er-DEXP<-X2-SNR) 

CALL  FNF11(-.5D0j  »5D0jX2»F11) 
AOR=E*Fll 

AOI=E*DSORT< .5D0*PI)*X 

CALL  FNF11 <~1 .5D0, .5D0»X2»F11 ) 

ANR=E*SNR#F 1 1 

ANI=SNR*<1 .5D0-X2  '*AOI 

FR=AOR+ANR 

FI=AOI+ANI 

BR=DMAX1 (DABS ( AOR) ? DABS <FR> ) 
BI=DMAX 1 (DABS(AOI ) »  DABS ( F I ) ) 

T= . 5D0+X2 

SNR2=SNR**2 
DO  100  N=2 1 200 
F0=N**2 

F1=SNR*(N+N-T)./F0 

F2  =  SNR2*<N-.5D0>/<  <N-1 )*F0> 

R=F1*ANR-F2*A0R 

V=F1*ANI-F2*A0I 

A0R~ANR 

AOI=ANI 

ANR  =  R 

ANI=V 

FR=FR+R 

FI*FI+V 

BR=DMAX1 ( BR » DABS  <  FR ) > 

BI=DMAX1 (BI»DABS(FI)) 

IF  <  DABS(U) .LE.5.D-19*DABS<FI ) 
GO  TO  200 
CONTINUE 
PRINT  101 

FORMAT <2X> '200  TERMS  IN  RICE') 
DR=18 . -DL0G10  <  DABS  <  BR/FR  > ) 
DI»18.-DL0G10< DABS (BI/FI) ) 
RETURN 


.AND.  DABS(R) . LE . 5 . D-19*DABS < FR >  > 


no  on  on  on  non  no 
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SUBROUTINE  CMPDVA 
PARAMETER  MF=2**10 
PARAMETER  F‘BDNUM=18 
DOUBLE  PRECISION  AL > AD * ABS * BSA * AM * ALPHA » ALFA » X ( 0 J 1 023 ) *Y (01 1023) 
PARAMETER  <NUMFIL=30) 

CHARACTER**  F ILES < NUMFIL > 

COMMON  /FILEC/FILES 
CHARACTER**  PBDNAM 
EQUIVALENCE 

1  <PBDNAM»FILES(18)  ) 

PARAMETER  <NUMPAR=200> 

COMMON  /PARAMC/PARAMS(NUMPAR) 

EQUIVALENCE 

1  (SMPMINf PARAMSt 187) > * 

1  (SNMIN»PARAMS< 184) ) *  <  SNMAX *  PAR AMS (185))*  ( SNDEL  »PARAMS( 186) ) 
COMMQN/PDVPF/AL*  AD » ABS 
DOUBLE  PRECISION  PI 
DATA  PI/3. 14159265358979324D0/ 


OPEN  THE  FILE 

CALL  OPNFIL  ( F’BDNUM  *  PBDNAM  > 
COMPUTE  THE  NUMBER  OF  SNR  CURVES 
NSN=<SNMAX-SNMIN)/SNDEL  +  i 


STORE  HEADER  INFO 

WR I TE ( PBDNUM )  SMPMIN  *  SNMIN  * SNMAX » SNDEL  *  NSN 
AM  -  SMPMIN 

AL  ••=  DMI N1 <  9 . DO  *  17 . DO/DSQRT  <  AM) ) 

AD  =  . 12D0/DSQRT ( AM ) 

BSA  ”  ~  DSQRT  <  P 1/2 . DO ) *AM  +  6 . DO*DSQRT ( AM ) 
ABS  =  DMIN1 (O.DO* BSA) 

COMPUTE  SNR  VS  PFA 
ALFARO. DO 

CALI  DIST  (AMjALFA*  MF »  X  »  Y ) 

STORE  THE  SNR  VS  PD 
WR II E ( PBDNUM )  <Y(I)»I =0*512) 

DO  1000  ISN=1*NSN 
SNR=  SNMIN  +  SNDEL*  <  I SN--.1  ) 

ALPHA  »  SNR 

CALL  DIST<AM*ALPHA*MF*X*Y> 

STORE  THE  SNR  VS  PD 
WR I TE  <  PBDNUM  >  ( Y ( I >  » 1=0*  51 2) 

1000  CONTINUE 

2000  CONTINUE 


RETURN 

END 
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SUBROUTINE  PRTSNR 

IMPLICIT  DOUBLE  PRECISION  <A--H,0-Z) 
DIMENSION  PFA(IO) fPD(lO) , V ( 1 4 > 8 >  *  SNR ( 1 4 , 8 ) 
DATA  PI/3.14159265358979324D0/ 

REAL  SMPMIN>SMPMAX*SMPFCT>PARAMS 
PARAMETER  NUMPAR=200 
COMMON/P ARAMC/P ARAMS  <NUMPAR) 

EQUIVALENCE 

1  ( SMPMIN »P ARAMS (187) ) >  < SMPMAX > PARAMS < 1 88 > )  • 

COMMON/PDPF /NPD»NPFAfPD»PFA 


MMAX  =  AL0G10(SMPMAX/SMPMIN>./AL0G10<SMPFCT>  + 


F 1 =DSQRT < ,5D0*PI ) 

F2=DSQRT  <2.D0-.5D0*PI) 

DO  1000  IM  =  1  *  MMAX 
AM=SMPMIN*SMPFCT**< IM-1 > 

AL  «  DMIN1 <9,B0» 1 7 . DO/DSQRT < AM > > 

AD  -  . 12D0/DSQRT<AM> 

BSA  =  -DSQRT<PI/2.nO)*AM  +  6»D0#DSQRT(AM) 
ABS  =  DMIN1(0.D0*  BSA ) 

AMU=F1*AM 

SIG=F2*DSQRT(AM) 

DO  900  IPF  =  1 » NPFA 
PF  =  PFA  < I PF ) 

IF (  AM  *  GT  *  1 ,  DO  )  GO  TO  250 
VN=DSQRT<-2.*DL0G<PF>  > 

GO  TO  750 

250  CALL  FNIPHI <  PF » YF ) 

V1=AMU-SIG*YF+ABS 

IF  (  I PF  *  GT  » 1  )  V1=DMAX1  < VI  * VN) 

V2=V1 + » 5D0 

IF (  Vl.NE.VN  )  GO  TO  300 
PI  =PN 
GO  TO  325 

300  CALL  FNPF<Vl»AM»AL»ADfABS»Pl) 

325  CALL  FNPF (V2»AM»AL*AD»ABS>P2) 

IF (  DABS (Pl-PF) »LT  *DABS<P2-PF)  )  GO  TO  350 
V0=V1 
P0  =  P1 
VN=V2 
PN  =  P2 
GO  TO  400 
350  V0=V2 

P0==P2 
VN=V1 
PN=P1 

400  CALL  FNIPHI (POtYO) 

GO  TO  550 

500  CALL  FNPF<VN*AM»AL»AD>  ABS  >  PN ) 

550  CALL  FNIPHI < PN » YN ) 

IF <  DABS<PN-PF) ,LE.1D-9*PF  )  GO  TO  750 
T=  <  VO#  <  YN-YF ) +VN* (YF-YO) )/<YN-YO) 


( SMPFCT  t PARAME 


1 
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VO=VN 
YO  =  YN 
VN  =  T 

GO  TO  500 
U<IM*IPF)=VN 
CONTINUE 
CONTINUE 


DO  4000  IPD-l » NPD 


CALL  FNIPHI < PD ( IPD) »  YD) 

DO  3000  IM--1»MMAX 
AM=SMPMIN*SMPFCT**< IM-i> 

AL  =  DMIN1<9*D0»17*D0/DSQRT<AM)) 

AD  =  . 12D0/DSGRT (AM) 

BSA  =  -DSQRT (PI/2  » DO ) #AH  +  6  * DO*DSQRT ( AM ) 

ADS  =  DMIN1 <0»D0»BSA> 

DO  2900  IPF=1 »NPFA 
PF=PFA< IFF) 

CALL  RDC  <  AM  r PF »PD(IPD)»A1> 

A2=A1*1 »01D0 
UV  =  U ( IM » IPF ) 

CALL  FNPD<A1»VU»AH»AL»  AD  * ABS » PI > 

CALL  FNPD ( A2 »  VV »  AM » AL  » AD » ABS » P2) 

IF (  DABS<P1-PD(IPD) ) »LT . DABS (P2-PD ( IPD) )  )  GO  TO  2350 
A0  =  A1 
P0  =  P1 
AN  =  A2 
PN  =  P2 

GO  TO  2400 
AO  =  A? 

P0=P2 
AN^Al 
PN  =  P1 

CALL  FNIPHI ( PO » YO ) 

GO  TO  2550 

CALL  FNPD(AN»VU»AM*AL>AD»ABS»PN) 

CALL  FNIPHI (PN>YN) 

IF (  DABS<PN-PD<IPD> ) »LE»1D-6#PD( IPD)  >  GO  TO  2750 
T~< A0*< YN-YD)+AN*< YD-YO) )/< YN-YO) 

AO  =  AN 
Y0~  YN 
AN~T 

GO  TO  2500 

SNR<IM>IPF)~10.#DL0G10< » 5D0*AN#AN ) 

CONTINUE 

CONTINUE 
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DO  3200  IPF=1jNPFA 
PRINT  3001 
3001  FORMAT ( 2  < / ) ) 

PRINT  301 1 f  PD<IPD)>PFA<IPF) 

3011  FORMAT  <  2X » 'PD  ~ ' *F10*3»5X» '  PFA  ~',D10.3> 
DO  3100  IM"1 »MMAX 
M=SMPMIN*SMPFCT**< IM-1  ) 

PRINT  3021 »  M  t SNR  < IM » IPF ) 

3021  FORMAT <  2X  > '  M  =  ' » 15  *  5X  > 'SNR  =',F7.2> 

3100  CONTINUE 
3200  CONTINUE 

4000  CONTINUE 


RETURN 

END 


SUBROUTINE  CMPSVS 

IMPLICIT  DOUBLE  PRECISION  (A-H..O-Z) 

PARAMETER  MMAX=14 

PARAMETER  NUMFIL=30>  PBDNUM=18 

CHARACTER*6  FILES<  NIJMFIL ) 

COMMON/FILEC/FILES 
CHARACTER*6  PBDNAM 
EQUIVALENCE  < PBDNAM » F ILES ( 1 8  )  > 

DIMENSION  PFA(10>  »PD<10)s>V(14»8)  s>  ALPHA  ( 1 4  » 8  > 
DIMENSION  THRS ( 1 4  »  8 ) 

DATA  PI/3.14159265358979324D0/ 
COMMON/PDPF/NPD » NPFA  >  PD  j>  PF  A 


CALL  OPNFIL ( PBDNUM » PBDNAM ) 


B-20 
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900 

1000 


F 1 =DSQRT ( .5D0*PI ) 

F2=DSGRT(2,D0-.5D0*PI> 

DO  1000  IM-1 »  MMAX 
AM  =  2 . #* ( IM-1 ) 

AL  =  DMIN1 <9»D0» 17  »DO/DSQRT (AM) ) 

AD  =  *  1 2D0/DSQRT  <  AM ) 

BSA  =  -DSQRT<PI/2.D0)#AM  +  6.D0#DSQRT(AM) 
ABS  =  DMIN1 <OtDO>BSA) 

AMU=F1*AM 
SIG=F2*DSGRT< AM) 

DO  900  1PF=1»8 
PF=10.**<-DFL0TJ< IPF) ) 

IF  <  AM ♦ GT  »  1  *  DO  )  GO  TO  250 
VN=DSGRT(-2.*DL0G<PF>  > 

GO  TO  750 

CALL  FNIPHI <  PF  >  YF ) 

V1=AMU-SIG*YF+ABS 

IF  <  IPF  »GT  ► 1  )  V1=DMAX1 (VI » VN ) 

V2=V1 +  ► 5D0 

IF (  VI . NE » ON  )  GO  TO  300 
PI  ~PN 
GO  TO  325 

CALL  FNPF(V1 , AM , AL » AD » ABS , P 1 ) 

CALL  FNPF ( V2>  AM » AL  » AD » ABS  » P2 ) 

IF (  DABS <P1-PF) »LT  »DABS<P2-PF>  )  GO  TO  350 

V0  =  V1 

P0  =  P1 

VN  =  V2 

PN  =  P2 

GO  TO  400 

V0=V2 

P0  =  P2 

VN  =  V1 

PN  =  F'l 

CALL  FNIPHI <P0> Y0> 

GO  TO  550 

CALL  FNPF(VN»AM>AL»ADjABS»PN) 

CALL  FNIPHI ( PN » YN ) 

IF  (  DABS<PN--PF>  ,LE.1D-9*PF  )  GO  TO  750 
T=( V0*< YN-YF)+VN#< YF-YO) )/< YN-YO) 

VO=VN 
YO -YN 
VN  =  T 

BO  10  500 
V< IM* IPF)-VN 

THRS  < IM » IPF >~<VN-ABS)/AM 

CONTINUE 

CONTINUE 


K>  fO  M  U 
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2350 


2400 

500 

550 


750 

900 

3000 


4000 


WRITE  (PBDNUM)  NPD>  (PD(I).MjIO) 

DO  4000  IPD=1jNFD 


CALL  FNIPHI <PD(1 PD) *  YD) 

DO  3000  IM-l  t  MMAX 
AM-2 ► DO*# ( IM- 1 ) 

AL.  =  DMIN1  ( 9  *  DO  >  17  *  D0/D3QRT  ( AM )  ) 

AD  =  . 12D0/DSQRT( AM) 

BSA  =  -DSQRT (PI/2*D0)#AM  +  6 . DOtDBQRT  <  AM ) 

ABS  -  DMIN1 <  0 . DO  »  BSA ) 

DO  2900  IPF~1 *  8 

PF=1 O.DO**(-DFLOTJ( IPF) ) 

CALL  RDC(AM>PF,PD(  IPD)  >  Al.  > 

A2=A1*1 >01D0 
VV  =  V ( IM  > IPF ) 

CALL  FNPD<  AJ  ,  VU ,  AM , AL  »  AD »  ABS  ,  P.t  ) 

CALL  FNPD ( A2  >  MU  *  AM  j AL  *  AD  * ABS  » P2 ) 

IF  <  DABS  ( PI  -  PD  <  IPD  )  )  .  LT  ♦  DABS  ( F'2~F’D  ( IF'D  )  )  )  GO  TO 

A0-A1 

P0=-P1 

AN~A2 

PN-P2 

GO  TO  2400 
A0--A2 
P0~F'2 
AN  =  A1 
PN  =  P1 

CALL  FNIPHI ( PO  >  YO ) 

GO  TO  2550 

CALL  FNFD< AN,VV> AM » AL » AD > ABS * PN > 

CALL  FNIPHI (PN>YN) 

IF (  DABS <PN-PD< IPD) ) . LE , 1 D-6*PD ( IPD )  )  GO  TO  2750 
T=< A0*< YN-YD)+AN*<YD-YO> )/( YN-YO) 

AO-AN 
YO^YN 
AN  =  T 

GO  TO  2500 
ALPHA ( IM » IPF ) =AN 
CONTINUE 
CONTINUE 


WRITE <PBDNUM)  < ( ALPHA < IM » IPF ) , IPF^l , 8 ) , IM=1 , MMAX ) 
CONTINUE 


RETURN 
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